#define DELTA_POS 0.00001
#include <algorithm>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <limits>
#include <boost\lexical_cast.hpp>
#include "wingbox.h"

wingbox::wingbox (int number_stringers, const stringer& _s, const box& _b, const rivet& _r, double rib1pos, double rib2pos, double loadbuckling, double loadfail, double safety)
	: b(_b),
	rib1(rib1pos / 1000.0),
	rib2(rib2pos / 1000.0),
	bucklingload(loadbuckling),
	failureload(loadfail),
	safetyfactorstringers(safety),
	safetyfactorrivets(safety),
	poprivets(_r)
{
	minimum_spacing = b.getWidth()  / (number_stringers-1);
	createTopStringers(_s);
	createBotStringers(_s);
	//cutStringers();

}
wingbox::wingbox (int number_stringers, const stringer& _s, const box& _b, const rivet& _r, double rib1pos, double rib2pos, double loadbuckling, double loadfail, double safety_stringers, double safety_rivets)
	: b(_b),
	rib1(rib1pos / 1000.0),
	rib2(rib2pos / 1000.0),
	bucklingload(loadbuckling),
	failureload(loadfail),
	safetyfactorstringers(safety_stringers),
	safetyfactorrivets(safety_rivets),
	poprivets(_r)
{
	minimum_spacing = b.getWidth()  / (number_stringers-1);
	createTopStringers(_s);
	createBotStringers(_s);
	cutStringers();

}
void wingbox::createTopStringers(const stringer& _s)
{
	char N('A');
	position p(1.5, /*_s.getWidth() / 2*/ 0, -b.getThickness(), -1, 0);
	for (; p.y < b.getWidth() - _s.getWidth(); p.y += minimum_spacing, ++p.index) {
		top_stringers.push_back(std::make_tuple(_s, p, std::vector<double>()));
		std::get<0>(top_stringers.back()).getName() = boost::lexical_cast<std::string>(N++);
	}
	p.y = b.getWidth() /*- _s.getWidth() / 2*/;
	top_stringers.push_back(std::make_tuple(_s, p, std::vector<double>()));
	std::get<0>(top_stringers.back()).getName() = boost::lexical_cast<std::string>(N++);
	createPermutations();
	cutTopStringers();

}
void wingbox::createBotStringers(const stringer& _s)
{
	position p(1.5,_s.getWidth() / 2,-b.getHeight() + b.getThickness());
	
	bottom_stringers.push_back(std::make_tuple(_s, p, std::vector<double>()));
	p.y = b.getHeight();
	++p.index;
	bottom_stringers.push_back(std::make_tuple(_s, p, std::vector<double>()));
	p.y = b.getWidth()- _s.getWidth() / 2;
	++p.index;
	bottom_stringers.push_back(std::make_tuple(_s, p, std::vector<double>()));
	cutBotStringers();
}
void wingbox::createPermutations() 
{
	unique_permutations.clear();
	auto orderfun( [] (const stringerlist::value_type& sdata1, const stringerlist::value_type& sdata2) {
		return std::get<1>(sdata1).y < std::get<1>(sdata2).y;
	});
	std::sort(top_stringers.begin()+1, top_stringers.end()-1, orderfun);
    do {
		removeorder_list l;
		for (auto b(top_stringers.begin()+1); b != top_stringers.end()-1; ++b) {
			l.push_back(std::get<1>(*b).index);
		}
		unique_permutations.push_back(l);
	} while(std::next_permutation(top_stringers.begin()+1, top_stringers.end()-1, orderfun));


	std::sort(top_stringers.begin()+1, top_stringers.end()-1, orderfun);
}

double wingbox::inertia() const
{
	double c(centroid());
	
	return b.inertia() +sqr(-b.centroid() - c) * b.area() 
		+ inertiaBotStringers() + sqr(centroidBotStringers() - c) * areaBotStringers()
		+ inertiaTopStringers() + sqr(centroidTopStringers() - c) * areaTopStringers();
}
double wingbox::centroid() const
{
	return (-b.centroid() * b.area() + centroidBotStringers() * areaBotStringers() + centroidTopStringers() * areaTopStringers()) /
		(area());
}
double wingbox::weight() const
{
	return b.weight() + weightTopStringers() + weightBotStringers() + weightRivets();
}
double wingbox::area() const
{
	return b.area() + areaBotStringers() + areaTopStringers() ;
}
double wingbox::inertia(double p) const
{
	double c(wingbox::centroid(p));
	
	return b.inertia() +sqr(-b.centroid() - c) * b.area() 
		+ inertiaBotStringers(p) + sqr(centroidBotStringers() - c) * areaBotStringers(p)
		+ inertiaTopStringers(p) + sqr(centroidTopStringers() - c) * areaTopStringers(p);
}
double wingbox::centroid(double p) const
{
	return (-b.centroid() * b.area() + centroidBotStringers() * areaBotStringers(p) + centroidTopStringers() * areaTopStringers(p)) /
		(area(p));
}
double wingbox::area(double p) const
{
	return b.area() + areaTopStringers(p) + areaBotStringers(p);
}


double wingbox::getKc(double pos) const
{
	return b.getKc();
}
double wingbox::getSweepFactor(double pos) const 
{
	if (pos  > rib1) {
		return 4.0/3.0;
	}
	return 1;
}

void wingbox::cutStringers()
{
	cutBotStringers();
	cutTopStringers();
}
void wingbox::cutTopStringers()
{
	cutStringers(top_stringers.begin(), top_stringers.end());
}
void wingbox::cutBotStringers()
{
	cutStringers(bottom_stringers.begin(), bottom_stringers.end());
}
void wingbox::cutStringers(stringer_iterator b, stringer_iterator e) 
{
	std::for_each(b, e, [this] (stringerlist::value_type& sdata) {
		cutStringer(sdata);
	});
}
void wingbox::cutStringer(stringerdata& sdata)
{
	position& p(std::get<1>(sdata));
	stringer& s(std::get<0>(sdata));
	double w = this->b.getWidth() - p.y;
	double d = w * std::tan(this->b.getSweep());
	double oldx = p.x;
	p.x = this->b.getLength() - d;
	double deltax = oldx - p.x;
	s.getLength() -= deltax;
}

double wingbox::centroidBotStringers() const
{
	return centroidStringers(bottom_stringers.front());
}
double wingbox::areaBotStringers() const
{
	return areaStringers(bottom_stringers.begin(), bottom_stringers.end());;
}
double wingbox::weightBotStringers() const
{
	return weightStringers(bottom_stringers.begin(), bottom_stringers.end());
}
double wingbox::centroidTopStringers() const
{
	return centroidStringers(top_stringers.front());
}
double wingbox::areaTopStringers() const
{
	return areaStringers(top_stringers.begin(), top_stringers.end());;
}
double wingbox::weightTopStringers() const
{
	return weightStringers(top_stringers.begin(), top_stringers.end());
}
double wingbox::areaStringers(stringer_const_iterator b, stringer_const_iterator e) const
{
	double stringer_area(0);
	std::for_each(b, e, [&stringer_area] (const stringerlist::value_type& sdata) {
		stringer_area += std::get<0>(sdata).area();
	});
	return stringer_area;
}
double wingbox::centroidStringers(const stringerlist::value_type& sdata) const
{
	const stringer& s(std::get<0>(sdata));
	const position& p(std::get<1>(sdata));
	return p.orientation * s.centroid() + p.z;
}
double wingbox::weightStringers(stringer_const_iterator b, stringer_const_iterator e) const
{
	double stringer_weight(0);
	std::for_each(b, e, [&stringer_weight] (const stringerlist::value_type& sdata) {
		stringer_weight += std::get<0>(sdata).weight();
	});
	return stringer_weight;
}
double wingbox::inertiaBotStringers() const
{
	return inertiaStringers(bottom_stringers.begin(), bottom_stringers.end());
}
double wingbox::inertiaTopStringers() const
{
	return inertiaStringers(top_stringers.begin(), top_stringers.end());
}
double wingbox::inertiaStringers(stringer_const_iterator b, stringer_const_iterator e) const
{
	double inertia(0);
	std::for_each(b, e, [&inertia] (const stringerlist::value_type& sdata) {
		inertia += std::get<0>(sdata).inertia();
	});
	return inertia;
}

double wingbox::areaBotStringers(double pos) const
{
	return areaStringers(bottom_stringers.begin(), bottom_stringers.end(), pos);
}
double wingbox::areaTopStringers(double pos) const
{
	return areaStringers(top_stringers.begin(), top_stringers.end(), pos);
}
double wingbox::areaStringers(stringer_const_iterator b, stringer_const_iterator e, double pos) const
{
	double stringer_area(0);
	std::for_each(b, e, [&stringer_area, &pos] (const stringerlist::value_type& sdata) {
		if (wingbox::stringerAtPos(sdata, pos)) {
			stringer_area += std::get<0>(sdata).area();
		}
	});
	return stringer_area;
}

double wingbox::inertiaBotStringers(double pos) const
{
	return inertiaStringers(bottom_stringers.begin(), bottom_stringers.end());
}
double wingbox::inertiaTopStringers(double pos) const
{
	return inertiaStringers(top_stringers.begin(), top_stringers.end());
}
double wingbox::inertiaStringers(stringer_const_iterator b, stringer_const_iterator e, double pos) const
{
	double inertia(0);
	std::for_each(b, e, [&inertia, &pos] (const stringerlist::value_type& sdata) {
		if (wingbox::stringerAtPos(sdata, pos)) {
			inertia += std::get<0>(sdata).inertia();
		}
	});
	return inertia;
}

bool wingbox::stringerAtPos(const stringerlist::value_type& sdata, double pos)
{
	double endpos( std::get<1>(sdata).x - std::get<0>(sdata).getLength() );
	return endpos + DELTA_POS < pos;
}

bool wingbox::setStringerSpacings(const std::vector<double>& spacings)
{
	std::vector<double> oldspaces;
	position prevpos = (std::get<1>(top_stringers[0]));
	std::for_each(top_stringers.begin() + 1, top_stringers.end(), [this, &oldspaces, &prevpos] (const stringerdata& sdata) {
		oldspaces.push_back(abs(std::get<1>(sdata).y - prevpos.y ));
		prevpos = std::get<1>(sdata);
	});


	int i(0);
	prevpos = (std::get<1>(top_stringers[0]));
	std::for_each(spacings.begin(), spacings.end(), [this, &prevpos, &i, &oldspaces] (double v) {
		double space;
		double maximumspacing = this->maximum_spacing(prevpos.x).second;
		if (v == -1) {
			space = maximumspacing;
		} else if (v == 0) {
			space = oldspaces[i];
		}else {
			space = std::min(v, maximumspacing);
		}
		++i;
		std::get<1>(this->top_stringers[i]).y = prevpos.y + space;
		cutStringer(top_stringers[i]);
		prevpos = std::get<1>(this->top_stringers[i]);
	});

	this->findBestRemovalOrder();
	return this->maximum_spacing(prevpos.x).second >= abs(prevpos.y - std::get<1>(top_stringers[++i]).y);
}
std::pair<int, double> wingbox::maximum_spacing(double pos) const
{
	double b_s(maximum_spacing_sheet_buckling(pos));
	return std::make_pair(1, b_s);
	//double b_c(maximum_spacing_column_buckling(pos));
	//if (b_s < b_c) {
	//	return std::make_pair(1, b_s);
	//} else {
	//	return std::make_pair(2, b_c);
	//}
}
double wingbox::maximum_spacing_column_buckling(double pos) const {
	return 0;
}

double wingbox::maximum_spacing_sheet_buckling(double pos) const
{
	return std::sqrt(  abs((  getKc(pos) * b.getMaterial().youngsmodulus * sqr(b.getThickness()) * this->inertia(pos)) 
		/ ( getSweepFactor(pos) * safetyfactorstringers * bucklingload * pos * this->centroid(pos)) ) );
}
double wingbox::findSpacingPossible(double spacing, double moment_of_inertia, double centroid_, double kc) const
{
	double T_over_B_squared = sqr(b.getThickness() / spacing);
	double v = ( kc * b.getMaterial().youngsmodulus * T_over_B_squared * moment_of_inertia) / ( getSweepFactor(1.5) * safetyfactorstringers * bucklingload * abs(centroid_));
	if (v <= rib1) {
		return findSpacingPossibleNoSweep(spacing, moment_of_inertia, abs(centroid_), kc);
	}
	return v;
}
double wingbox::findSpacingPossibleNoSweep(double spacing, double moment_of_inertia, double centroid_, double kc) const
{
	double T_over_B_squared = sqr(b.getThickness() / spacing);
	double v = ( kc * b.getMaterial().youngsmodulus * T_over_B_squared * moment_of_inertia) / ( safetyfactorstringers * bucklingload * abs(centroid_));
	return v;
}
void wingbox::displaydata() const
{
	std::cout << "wingbox: " << "" <<std::endl;
	std::cout << "stringers at top: " << top_stringers.size() << std::endl;
	std::cout << "stringers at bot: " << bottom_stringers.size() << std::endl;
	std::cout << std::setprecision(3) << "minimum spacing possible at root: " << this->minimum_spacing *1000 << std::endl;
	std::cout << std::setprecision(3) << "maximum spacing possible at root (leading edge): " << maximum_spacing(b.getLength()).second *1000 <<std::endl;
	std::cout << std::setprecision(3) << "maximum spacing possible at root (trailing edge): " << maximum_spacing(b.getLength() - b.getWidth() * std::tan(b.getSweep())).second *1000 <<std::endl;
	std::cout << std::setprecision(0) << std::endl;
	std::cout << "centroid: " << centroid() << std::endl 
			<< "moment of inertia: " << inertia() <<std::endl
			<< "cross section area: " << area() << std::endl << std::endl;

	std::cout << "total mass: " << weight() << std::endl << std::endl;
	
	displaystringerdata();


}
void wingbox::displaystringerdata() const
{
	std::cout << "stringer pos & length: " <<std::endl;
	std::for_each(top_stringers.begin(), top_stringers.end(), [this] (const stringerlist::value_type& sdata) {
		const position& p(std::get<1>(sdata));
		const stringer& s(std::get<0>(sdata));
		std::cout << s.getName() << std::setprecision(3) << ": { "  << (1.5 - p.x) * 1000  << " , "  << p.y * 1000 << " } l: " 
			<< std::setprecision(4)  << s.getLength() * 1000 << std::setprecision(0) << std::endl;
	});


}
int wingbox::numberTopStringers() const
{
	return top_stringers.size();
}

void wingbox::findBestRemovalOrder()
{
	
	permutation_list::const_iterator best_permutation(unique_permutations.end());
	std::vector<double> origlengths;
	std::vector<double> bestlengths;
	std::for_each(top_stringers.begin()+1, top_stringers.end()-1, [&origlengths, this] (stringerlist::value_type& v) {
		std::get<0>(v).getLength() = std::get<1>(v).x;

		origlengths.push_back(std::get<0>(v).getLength());
	});
	auto maxweight = weight();
	auto e(unique_permutations.end());
	for (auto permutation(unique_permutations.begin()); permutation != e; ++permutation) {
		double w = testRemovalOrder(*(permutation));
		if (w < maxweight) {
			best_permutation = (permutation);
			bestlengths.clear();
			std::for_each(top_stringers.begin()+1, top_stringers.end()-1, [&bestlengths, this] (const stringerlist::value_type& v) {
				bestlengths.push_back(std::get<0>(v).getLength());
			});
			maxweight = w;
		}
		auto origlength_end(origlengths.end());
		int i(1);
		for (auto iter_origlength(origlengths.begin());  iter_origlength != origlength_end; ++iter_origlength, ++i) {
			std::get<0>(top_stringers[i]).getLength() = *iter_origlength;
		}
	}

	auto bestlengths_end(bestlengths.end());
	int i(1);
	std::cout << "best option: ";
	std::for_each(best_permutation->begin(), best_permutation->end(), [this] (const decltype(*best_permutation->begin()) v) {
		std::cout << std::get<0>(top_stringers[v]).getName();
	});
	std::cout << std::endl;
	for (auto iter_bestlengths(bestlengths.begin());  iter_bestlengths != bestlengths_end; ++iter_bestlengths, ++i) {
		std::get<0>(top_stringers[i]).getLength() = *iter_bestlengths +0.005;
	}
	
}
double wingbox::calcSpaceWithout(wingbox::stringer_const_iterator i, double pos) const 
{
	auto j(i - 1);
	while (!wingbox::stringerAtPos(*j, pos) && j != top_stringers.begin()) {
		--j;
	}
	++i;
	while (i != top_stringers.end() -1 && !wingbox::stringerAtPos(*i, pos) ) {
		++i;
	}
	return abs(std::get<1>(*j).y - std::get<1>(*i).y);
}
double wingbox::calcSpace(double xpos, double ypos)
{
	auto iter(top_stringers.begin());
	auto first(iter);
	do {
		do {
			first = iter;
			if (iter != top_stringers.end()) {
				return -1;
			}
			++iter;
		} while ( std::get<1>(*iter).x - std::get<0>(*iter).getLength() > xpos) ;
	} while (std::get<1>(*iter).y < ypos);
	
	do {
		if (iter != top_stringers.end()) {
			return -1;
		}
		++iter;
	} while ( std::get<1>(*iter).x - std::get<0>(*iter).getLength() > xpos) ;
	return std::get<1>(*iter).y - std::get<1>(*first).y;

}

double wingbox::testRemovalOrder(const removeorder_list& permutation) {
	auto maxweight = std::numeric_limits<double>::max();
	double pos(b.getLength());
	//std::cout << "Removal order: ";
	//std::string debugstr = "";
	//std::for_each(permutation.begin(), permutation.end(), [this, &pos, &debugstr] (int i) {
	//	debugstr += std::get<0>(this->top_stringers[i]).getName();
	//});

	std::for_each(permutation.begin(), permutation.end(), [this, &pos] (int i) {
		//simulate away:
		std::string name(std::get<0>(top_stringers[i]).getName());
		std::get<0>(this->top_stringers[i]).getLength() = 0;
		auto titer(top_stringers.begin());
		std::advance(titer,i);
		double b(calcSpaceWithout(titer, pos));
		double I(inertia(pos));
		double c(centroid(pos));
		double kc(getKc(pos));
		double x;
		if (pos <= rib1) {
			x = findSpacingPossibleNoSweep(b, I, c, kc);
		} else {
			x = findSpacingPossible(b, I, c, kc);
		}
		if (x > pos) {
			x = pos;
		} else {
			pos = x;
		}
		std::get<0>(this->top_stringers[i]).getLength() = abs(std::get<1>(top_stringers[i]).x - x);
	});
	//std::cout  << " : " << weight() << '\n';
	return weight();
}

double& wingbox::getSafetyStringers()
{
	return safetyfactorstringers;
}
const double& wingbox::getSafetyStringers() const
{
	return safetyfactorstringers;
}
double& wingbox::getSafetyRivets()
{
	return safetyfactorrivets;
}
const double& wingbox::getSafetyRivets() const
{
	return safetyfactorrivets;
}

void wingbox::changeMaterial(const material& m)
{
	b.getMaterial() = m;
	std::for_each(top_stringers.begin(), top_stringers.end(), [&m] (stringerdata& sdata) {
		std::get<0>(sdata).getMaterial() = m;
	});
	std::for_each(bottom_stringers.begin(), bottom_stringers.end(), [&m] (stringerdata& sdata) {
		std::get<0>(sdata).getMaterial() = m;
	});
}

void wingbox::changeBox(double w, double h, double t, double l, double s)
{
	b.getWidth() = w / 1000;
	b.getHeight() = h / 1000;
	b.getThickness() = t / 1000;
	b.getLength() = l / 1000;
	b.setSweep(s);
}
void wingbox::changeStringer(int n, double w, double t)
{
	

	std::for_each(bottom_stringers.begin(), bottom_stringers.end(), [&] (stringerdata& sdata) {
		std::get<0>(sdata).getWidth() = w / 1000;
		std::get<0>(sdata).getThickness() = t / 1000;
	});

	if (n == top_stringers.size()) {
		std::for_each(top_stringers.begin(), top_stringers.end(), [&] (stringerdata& sdata) {
			std::get<0>(sdata).getWidth() = w / 1000;
			std::get<0>(sdata).getThickness() = t / 1000;
		});
	} else {
		
		stringer s(w, t, b.getLength() * 1000, std::get<0>(top_stringers.front()).getMaterial());
		top_stringers.clear();
		minimum_spacing = b.getWidth() / (n - 1);
		createTopStringers(s);

	}
}
void wingbox::changeRibPositions(double rib1, double rib2)
{
	this->rib1 = rib1;
	this->rib2 = rib2;
}
void wingbox::changeSafetyFactors(double stringers, double rivets)
{
	safetyfactorstringers = stringers;
	safetyfactorrivets = rivets;
}
void wingbox::changeMaximumLoadings(double buckling_load, double ultimate_load)
{
	bucklingload = buckling_load;
	failureload = ultimate_load;
}
void wingbox::changeBucklingConstants(double Kc, double endfixity)
{
	b.getKc() = Kc;
	poprivets.getEndFixity() = endfixity;
}
void wingbox::changeClamps(double clamp_root, double clamp_tip)
{
	b.setRootClamp(clamp_root);
	b.setTipClamp(clamp_tip);
}

void wingbox::changeName(const std::string& newname)
{
	name = newname;
}